home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlaln2.f < prev    next >
Text File  |  1996-07-19  |  17KB  |  509 lines
  1.       SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
  2.      $                   LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
  3. *
  4. *  -- LAPACK auxiliary routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     October 31, 1992
  8. *
  9. *     .. Scalar Arguments ..
  10.       LOGICAL            LTRANS
  11.       INTEGER            INFO, LDA, LDB, LDX, NA, NW
  12.       DOUBLE PRECISION   CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
  13. *     ..
  14. *     .. Array Arguments ..
  15.       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), X( LDX, * )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  DLALN2 solves a system of the form  (ca A - w D ) X = s B
  22. *  or (ca A' - w D) X = s B   with possible scaling ("s") and
  23. *  perturbation of A.  (A' means A-transpose.)
  24. *
  25. *  A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
  26. *  real diagonal matrix, w is a real or complex value, and X and B are
  27. *  NA x 1 matrices -- real if w is real, complex if w is complex.  NA
  28. *  may be 1 or 2.
  29. *
  30. *  If w is complex, X and B are represented as NA x 2 matrices,
  31. *  the first column of each being the real part and the second
  32. *  being the imaginary part.
  33. *
  34. *  "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
  35. *  so chosen that X can be computed without overflow.  X is further
  36. *  scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
  37. *  than overflow.
  38. *
  39. *  If both singular values of (ca A - w D) are less than SMIN,
  40. *  SMIN*identity will be used instead of (ca A - w D).  If only one
  41. *  singular value is less than SMIN, one element of (ca A - w D) will be
  42. *  perturbed enough to make the smallest singular value roughly SMIN.
  43. *  If both singular values are at least SMIN, (ca A - w D) will not be
  44. *  perturbed.  In any case, the perturbation will be at most some small
  45. *  multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values
  46. *  are computed by infinity-norm approximations, and thus will only be
  47. *  correct to a factor of 2 or so.
  48. *
  49. *  Note: all input quantities are assumed to be smaller than overflow
  50. *  by a reasonable factor.  (See BIGNUM.)
  51. *
  52. *  Arguments
  53. *  ==========
  54. *
  55. *  LTRANS  (input) LOGICAL
  56. *          =.TRUE.:  A-transpose will be used.
  57. *          =.FALSE.: A will be used (not transposed.)
  58. *
  59. *  NA      (input) INTEGER
  60. *          The size of the matrix A.  It may (only) be 1 or 2.
  61. *
  62. *  NW      (input) INTEGER
  63. *          1 if "w" is real, 2 if "w" is complex.  It may only be 1
  64. *          or 2.
  65. *
  66. *  SMIN    (input) DOUBLE PRECISION
  67. *          The desired lower bound on the singular values of A.  This
  68. *          should be a safe distance away from underflow or overflow,
  69. *          say, between (underflow/machine precision) and  (machine
  70. *          precision * overflow ).  (See BIGNUM and ULP.)
  71. *
  72. *  CA      (input) DOUBLE PRECISION
  73. *          The coefficient c, which A is multiplied by.
  74. *
  75. *  A       (input) DOUBLE PRECISION array, dimension (LDA,NA)
  76. *          The NA x NA matrix A.
  77. *
  78. *  LDA     (input) INTEGER
  79. *          The leading dimension of A.  It must be at least NA.
  80. *
  81. *  D1      (input) DOUBLE PRECISION
  82. *          The 1,1 element in the diagonal matrix D.
  83. *
  84. *  D2      (input) DOUBLE PRECISION
  85. *          The 2,2 element in the diagonal matrix D.  Not used if NW=1.
  86. *
  87. *  B       (input) DOUBLE PRECISION array, dimension (LDB,NW)
  88. *          The NA x NW matrix B (right-hand side).  If NW=2 ("w" is
  89. *          complex), column 1 contains the real part of B and column 2
  90. *          contains the imaginary part.
  91. *
  92. *  LDB     (input) INTEGER
  93. *          The leading dimension of B.  It must be at least NA.
  94. *
  95. *  WR      (input) DOUBLE PRECISION
  96. *          The real part of the scalar "w".
  97. *
  98. *  WI      (input) DOUBLE PRECISION
  99. *          The imaginary part of the scalar "w".  Not used if NW=1.
  100. *
  101. *  X       (output) DOUBLE PRECISION array, dimension (LDX,NW)
  102. *          The NA x NW matrix X (unknowns), as computed by DLALN2.
  103. *          If NW=2 ("w" is complex), on exit, column 1 will contain
  104. *          the real part of X and column 2 will contain the imaginary
  105. *          part.
  106. *
  107. *  LDX     (input) INTEGER
  108. *          The leading dimension of X.  It must be at least NA.
  109. *
  110. *  SCALE   (output) DOUBLE PRECISION
  111. *          The scale factor that B must be multiplied by to insure
  112. *          that overflow does not occur when computing X.  Thus,
  113. *          (ca A - w D) X  will be SCALE*B, not B (ignoring
  114. *          perturbations of A.)  It will be at most 1.
  115. *
  116. *  XNORM   (output) DOUBLE PRECISION
  117. *          The infinity-norm of X, when X is regarded as an NA x NW
  118. *          real matrix.
  119. *
  120. *  INFO    (output) INTEGER
  121. *          An error flag.  It will be set to zero if no error occurs,
  122. *          a negative number if an argument is in error, or a positive
  123. *          number if  ca A - w D  had to be perturbed.
  124. *          The possible values are:
  125. *          = 0: No error occurred, and (ca A - w D) did not have to be
  126. *                 perturbed.
  127. *          = 1: (ca A - w D) had to be perturbed to make its smallest
  128. *               (or only) singular value greater than SMIN.
  129. *          NOTE: In the interests of speed, this routine does not
  130. *                check the inputs for errors.
  131. *
  132. * =====================================================================
  133. *
  134. *     .. Parameters ..
  135.       DOUBLE PRECISION   ZERO, ONE
  136.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  137.       DOUBLE PRECISION   TWO
  138.       PARAMETER          ( TWO = 2.0D0 )
  139. *     ..
  140. *     .. Local Scalars ..
  141.       INTEGER            ICMAX, J
  142.       DOUBLE PRECISION   BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21,
  143.      $                   CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21,
  144.      $                   LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R,
  145.      $                   UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S,
  146.      $                   UR22, XI1, XI2, XR1, XR2
  147. *     ..
  148. *     .. Local Arrays ..
  149.       LOGICAL            RSWAP( 4 ), ZSWAP( 4 )
  150.       INTEGER            IPIVOT( 4, 4 )
  151.       DOUBLE PRECISION   CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
  152. *     ..
  153. *     .. External Functions ..
  154.       DOUBLE PRECISION   DLAMCH
  155.       EXTERNAL           DLAMCH
  156. *     ..
  157. *     .. External Subroutines ..
  158.       EXTERNAL           DLADIV
  159. *     ..
  160. *     .. Intrinsic Functions ..
  161.       INTRINSIC          ABS, MAX
  162. *     ..
  163. *     .. Equivalences ..
  164.       EQUIVALENCE        ( CI( 1, 1 ), CIV( 1 ) ),
  165.      $                   ( CR( 1, 1 ), CRV( 1 ) )
  166. *     ..
  167. *     .. Data statements ..
  168.       DATA               ZSWAP / .FALSE., .FALSE., .TRUE., .TRUE. /
  169.       DATA               RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. /
  170.       DATA               IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4,
  171.      $                   3, 2, 1 /
  172. *     ..
  173. *     .. Executable Statements ..
  174. *
  175. *     Compute BIGNUM
  176. *
  177.       SMLNUM = TWO*DLAMCH( 'Safe minimum' )
  178.       BIGNUM = ONE / SMLNUM
  179.       SMINI = MAX( SMIN, SMLNUM )
  180. *
  181. *     Don't check for input errors
  182. *
  183.       INFO = 0
  184. *
  185. *     Standard Initializations
  186. *
  187.       SCALE = ONE
  188. *
  189.       IF( NA.EQ.1 ) THEN
  190. *
  191. *        1 x 1  (i.e., scalar) system   C X = B
  192. *
  193.          IF( NW.EQ.1 ) THEN
  194. *
  195. *           Real 1x1 system.
  196. *
  197. *           C = ca A - w D
  198. *
  199.             CSR = CA*A( 1, 1 ) - WR*D1
  200.             CNORM = ABS( CSR )
  201. *
  202. *           If | C | < SMINI, use C = SMINI
  203. *
  204.             IF( CNORM.LT.SMINI ) THEN
  205.                CSR = SMINI
  206.                CNORM = SMINI
  207.                INFO = 1
  208.             END IF
  209. *
  210. *           Check scaling for  X = B / C
  211. *
  212.             BNORM = ABS( B( 1, 1 ) )
  213.             IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
  214.                IF( BNORM.GT.BIGNUM*CNORM )
  215.      $            SCALE = ONE / BNORM
  216.             END IF
  217. *
  218. *           Compute X
  219. *
  220.             X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR
  221.             XNORM = ABS( X( 1, 1 ) )
  222.          ELSE
  223. *
  224. *           Complex 1x1 system (w is complex)
  225. *
  226. *           C = ca A - w D
  227. *
  228.             CSR = CA*A( 1, 1 ) - WR*D1
  229.             CSI = -WI*D1
  230.             CNORM = ABS( CSR ) + ABS( CSI )
  231. *
  232. *           If | C | < SMINI, use C = SMINI
  233. *
  234.             IF( CNORM.LT.SMINI ) THEN
  235.                CSR = SMINI
  236.                CSI = ZERO
  237.                CNORM = SMINI
  238.                INFO = 1
  239.             END IF
  240. *
  241. *           Check scaling for  X = B / C
  242. *
  243.             BNORM = ABS( B( 1, 1 ) ) + ABS( B( 1, 2 ) )
  244.             IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
  245.                IF( BNORM.GT.BIGNUM*CNORM )
  246.      $            SCALE = ONE / BNORM
  247.             END IF
  248. *
  249. *           Compute X
  250. *
  251.             CALL DLADIV( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI,
  252.      $                   X( 1, 1 ), X( 1, 2 ) )
  253.             XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
  254.          END IF
  255. *
  256.       ELSE
  257. *
  258. *        2x2 System
  259. *
  260. *        Compute the real part of  C = ca A - w D  (or  ca A' - w D )
  261. *
  262.          CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1
  263.          CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2
  264.          IF( LTRANS ) THEN
  265.             CR( 1, 2 ) = CA*A( 2, 1 )
  266.             CR( 2, 1 ) = CA*A( 1, 2 )
  267.          ELSE
  268.             CR( 2, 1 ) = CA*A( 2, 1 )
  269.             CR( 1, 2 ) = CA*A( 1, 2 )
  270.          END IF
  271. *
  272.          IF( NW.EQ.1 ) THEN
  273. *
  274. *           Real 2x2 system  (w is real)
  275. *
  276. *           Find the largest element in C
  277. *
  278.             CMAX = ZERO
  279.             ICMAX = 0
  280. *
  281.             DO 10 J = 1, 4
  282.                IF( ABS( CRV( J ) ).GT.CMAX ) THEN
  283.                   CMAX = ABS( CRV( J ) )
  284.                   ICMAX = J
  285.                END IF
  286.    10       CONTINUE
  287. *
  288. *           If norm(C) < SMINI, use SMINI*identity.
  289. *
  290.             IF( CMAX.LT.SMINI ) THEN
  291.                BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 2, 1 ) ) )
  292.                IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
  293.                   IF( BNORM.GT.BIGNUM*SMINI )
  294.      $               SCALE = ONE / BNORM
  295.                END IF
  296.                TEMP = SCALE / SMINI
  297.                X( 1, 1 ) = TEMP*B( 1, 1 )
  298.                X( 2, 1 ) = TEMP*B( 2, 1 )
  299.                XNORM = TEMP*BNORM
  300.                INFO = 1
  301.                RETURN
  302.             END IF
  303. *
  304. *           Gaussian elimination with complete pivoting.
  305. *
  306.             UR11 = CRV( ICMAX )
  307.             CR21 = CRV( IPIVOT( 2, ICMAX ) )
  308.             UR12 = CRV( IPIVOT( 3, ICMAX ) )
  309.             CR22 = CRV( IPIVOT( 4, ICMAX ) )
  310.             UR11R = ONE / UR11
  311.             LR21 = UR11R*CR21
  312.             UR22 = CR22 - UR12*LR21
  313. *
  314. *           If smaller pivot < SMINI, use SMINI
  315. *
  316.             IF( ABS( UR22 ).LT.SMINI ) THEN
  317.                UR22 = SMINI
  318.                INFO = 1
  319.             END IF
  320.             IF( RSWAP( ICMAX ) ) THEN
  321.                BR1 = B( 2, 1 )
  322.                BR2 = B( 1, 1 )
  323.             ELSE
  324.                BR1 = B( 1, 1 )
  325.                BR2 = B( 2, 1 )
  326.             END IF
  327.             BR2 = BR2 - LR21*BR1
  328.             BBND = MAX( ABS( BR1*( UR22*UR11R ) ), ABS( BR2 ) )
  329.             IF( BBND.GT.ONE .AND. ABS( UR22 ).LT.ONE ) THEN
  330.                IF( BBND.GE.BIGNUM*ABS( UR22 ) )
  331.      $            SCALE = ONE / BBND
  332.             END IF
  333. *
  334.             XR2 = ( BR2*SCALE ) / UR22
  335.             XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 )
  336.             IF( ZSWAP( ICMAX ) ) THEN
  337.                X( 1, 1 ) = XR2
  338.                X( 2, 1 ) = XR1
  339.             ELSE
  340.                X( 1, 1 ) = XR1
  341.                X( 2, 1 ) = XR2
  342.             END IF
  343.             XNORM = MAX( ABS( XR1 ), ABS( XR2 ) )
  344. *
  345. *           Further scaling if  norm(A) norm(X) > overflow
  346. *
  347.             IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN
  348.                IF( XNORM.GT.BIGNUM / CMAX ) THEN
  349.                   TEMP = CMAX / BIGNUM
  350.                   X( 1, 1 ) = TEMP*X( 1, 1 )
  351.                   X( 2, 1 ) = TEMP*X( 2, 1 )
  352.                   XNORM = TEMP*XNORM
  353.                   SCALE = TEMP*SCALE
  354.                END IF
  355.             END IF
  356.          ELSE
  357. *
  358. *           Complex 2x2 system  (w is complex)
  359. *
  360. *           Find the largest element in C
  361. *
  362.             CI( 1, 1 ) = -WI*D1
  363.             CI( 2, 1 ) = ZERO
  364.             CI( 1, 2 ) = ZERO
  365.             CI( 2, 2 ) = -WI*D2
  366.             CMAX = ZERO
  367.             ICMAX = 0
  368. *
  369.             DO 20 J = 1, 4
  370.                IF( ABS( CRV( J ) )+ABS( CIV( J ) ).GT.CMAX ) THEN
  371.                   CMAX = ABS( CRV( J ) ) + ABS( CIV( J ) )
  372.                   ICMAX = J
  373.                END IF
  374.    20       CONTINUE
  375. *
  376. *           If norm(C) < SMINI, use SMINI*identity.
  377. *
  378.             IF( CMAX.LT.SMINI ) THEN
  379.                BNORM = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ),
  380.      $                 ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) )
  381.                IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
  382.                   IF( BNORM.GT.BIGNUM*SMINI )
  383.      $               SCALE = ONE / BNORM
  384.                END IF
  385.                TEMP = SCALE / SMINI
  386.                X( 1, 1 ) = TEMP*B( 1, 1 )
  387.                X( 2, 1 ) = TEMP*B( 2, 1 )
  388.                X( 1, 2 ) = TEMP*B( 1, 2 )
  389.                X( 2, 2 ) = TEMP*B( 2, 2 )
  390.                XNORM = TEMP*BNORM
  391.                INFO = 1
  392.                RETURN
  393.             END IF
  394. *
  395. *           Gaussian elimination with complete pivoting.
  396. *
  397.             UR11 = CRV( ICMAX )
  398.             UI11 = CIV( ICMAX )
  399.             CR21 = CRV( IPIVOT( 2, ICMAX ) )
  400.             CI21 = CIV( IPIVOT( 2, ICMAX ) )
  401.             UR12 = CRV( IPIVOT( 3, ICMAX ) )
  402.             UI12 = CIV( IPIVOT( 3, ICMAX ) )
  403.             CR22 = CRV( IPIVOT( 4, ICMAX ) )
  404.             CI22 = CIV( IPIVOT( 4, ICMAX ) )
  405.             IF( ICMAX.EQ.1 .OR. ICMAX.EQ.4 ) THEN
  406. *
  407. *              Code when off-diagonals of pivoted C are real
  408. *
  409.                IF( ABS( UR11 ).GT.ABS( UI11 ) ) THEN
  410.                   TEMP = UI11 / UR11
  411.                   UR11R = ONE / ( UR11*( ONE+TEMP**2 ) )
  412.                   UI11R = -TEMP*UR11R
  413.                ELSE
  414.                   TEMP = UR11 / UI11
  415.                   UI11R = -ONE / ( UI11*( ONE+TEMP**2 ) )
  416.                   UR11R = -TEMP*UI11R
  417.                END IF
  418.                LR21 = CR21*UR11R
  419.                LI21 = CR21*UI11R
  420.                UR12S = UR12*UR11R
  421.                UI12S = UR12*UI11R
  422.                UR22 = CR22 - UR12*LR21
  423.                UI22 = CI22 - UR12*LI21
  424.             ELSE
  425. *
  426. *              Code when diagonals of pivoted C are real
  427. *
  428.                UR11R = ONE / UR11
  429.                UI11R = ZERO
  430.                LR21 = CR21*UR11R
  431.                LI21 = CI21*UR11R
  432.                UR12S = UR12*UR11R
  433.                UI12S = UI12*UR11R
  434.                UR22 = CR22 - UR12*LR21 + UI12*LI21
  435.                UI22 = -UR12*LI21 - UI12*LR21
  436.             END IF
  437.             U22ABS = ABS( UR22 ) + ABS( UI22 )
  438. *
  439. *           If smaller pivot < SMINI, use SMINI
  440. *
  441.             IF( U22ABS.LT.SMINI ) THEN
  442.                UR22 = SMINI
  443.                UI22 = ZERO
  444.                INFO = 1
  445.             END IF
  446.             IF( RSWAP( ICMAX ) ) THEN
  447.                BR2 = B( 1, 1 )
  448.                BR1 = B( 2, 1 )
  449.                BI2 = B( 1, 2 )
  450.                BI1 = B( 2, 2 )
  451.             ELSE
  452.                BR1 = B( 1, 1 )
  453.                BR2 = B( 2, 1 )
  454.                BI1 = B( 1, 2 )
  455.                BI2 = B( 2, 2 )
  456.             END IF
  457.             BR2 = BR2 - LR21*BR1 + LI21*BI1
  458.             BI2 = BI2 - LI21*BR1 - LR21*BI1
  459.             BBND = MAX( ( ABS( BR1 )+ABS( BI1 ) )*
  460.      $             ( U22ABS*( ABS( UR11R )+ABS( UI11R ) ) ),
  461.      $             ABS( BR2 )+ABS( BI2 ) )
  462.             IF( BBND.GT.ONE .AND. U22ABS.LT.ONE ) THEN
  463.                IF( BBND.GE.BIGNUM*U22ABS ) THEN
  464.                   SCALE = ONE / BBND
  465.                   BR1 = SCALE*BR1
  466.                   BI1 = SCALE*BI1
  467.                   BR2 = SCALE*BR2
  468.                   BI2 = SCALE*BI2
  469.                END IF
  470.             END IF
  471. *
  472.             CALL DLADIV( BR2, BI2, UR22, UI22, XR2, XI2 )
  473.             XR1 = UR11R*BR1 - UI11R*BI1 - UR12S*XR2 + UI12S*XI2
  474.             XI1 = UI11R*BR1 + UR11R*BI1 - UI12S*XR2 - UR12S*XI2
  475.             IF( ZSWAP( ICMAX ) ) THEN
  476.                X( 1, 1 ) = XR2
  477.                X( 2, 1 ) = XR1
  478.                X( 1, 2 ) = XI2
  479.                X( 2, 2 ) = XI1
  480.             ELSE
  481.                X( 1, 1 ) = XR1
  482.                X( 2, 1 ) = XR2
  483.                X( 1, 2 ) = XI1
  484.                X( 2, 2 ) = XI2
  485.             END IF
  486.             XNORM = MAX( ABS( XR1 )+ABS( XI1 ), ABS( XR2 )+ABS( XI2 ) )
  487. *
  488. *           Further scaling if  norm(A) norm(X) > overflow
  489. *
  490.             IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN
  491.                IF( XNORM.GT.BIGNUM / CMAX ) THEN
  492.                   TEMP = CMAX / BIGNUM
  493.                   X( 1, 1 ) = TEMP*X( 1, 1 )
  494.                   X( 2, 1 ) = TEMP*X( 2, 1 )
  495.                   X( 1, 2 ) = TEMP*X( 1, 2 )
  496.                   X( 2, 2 ) = TEMP*X( 2, 2 )
  497.                   XNORM = TEMP*XNORM
  498.                   SCALE = TEMP*SCALE
  499.                END IF
  500.             END IF
  501.          END IF
  502.       END IF
  503. *
  504.       RETURN
  505. *
  506. *     End of DLALN2
  507. *
  508.       END
  509.